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Abstract — We  present  a  computational  geometry  method  for 
the  problem  of  triangulation  in  the  plane  using  measurements 
of  distance-differences.  Compared  to  existing  solutions  to  this 
well-studied  problem,  this  method  is:  (a)  computationally  more 
efficient  and  adaptive  in  that  its  precision  can  be  controlled  as 
a  function  of  the  number  of  computational  operations,  making 
it  suitable  to  low  power  devices,  and  (b)  robust  with  respect  to 
measurement  and  computational  errors,  and  is  not  susceptible 
to  numerical  instabilities  typical  of  existing  linear  algebraic  or 
quadratic  methods.  This  method  employs  a  binary  search  on  a 
distance-difference  curve  in  the  plane  using  a  second  distance- 
difference  as  the  objective  function.  We  establish  the  unimodality 
of  the  directional  derivative  of  the  objective  function  within  each 
of  a  small  number  of  suitably  decomposed  regions  of  the  plane 
to  support  the  binary  search.  The  computational  complexity 
of  this  method  is  0(log2  1/7),  where  the  computed  solution  is 
guaranteed  to  be  within  a  7-precision  region  centered  at  the 
actual  solution.  We  present  simulation  results  to  compare  this 
method  with  existing  DTOA  triangulation  methods. 

Keywords:  Triangulation,  difference  in  time  of  arrival, 
computational  geometry,  computational  complexity. 

I.  Introduction 

The  problem  of  computing  the  location  of  an  object  from 
measurements  of  distance-differences  from  three  known  lo¬ 
cations  is  well-studied  (for  at  least  three  decades)  under  the 
title  of  Difference  of  Time-of-Arrival  (DTOA)  localization. 
This  problem  arises  in  a  number  of  established  areas  such  as 
tracking  in  aerospace  systems  [1],  [2].  Recently,  it  has  received 
renewed  attention  due  to  the  increasing  proliferation  of  wire¬ 
less  sensor  networks  [3],  [4]  and  embedded  networked  systems 
[5].  In  several  of  these  applications,  the  wireless  nodes  are  lim¬ 
ited  in  power  and  yet  the  localization  computations  may  have 
to  be  repeated  quite  frequently.  Consequently,  it  has  become 
important  to  trade-off  the  number  and  type  of  computations 
needed  for  localization  to  save  power  by  gracefully  degrading 
the  quality  of  solution.  In  addition,  the  computational  precision 
of  arithmetic  operations  may  be  limited  in  some  sensor  nodes, 
but  its  impact  on  the  precision  of  localization  is  not  well 
understood.  These  factors  motivate  a  closer  examination  of 
the  computational  aspects  of  DTOA  triangulation  methods; 
however,  our  results  could  be  of  more  general  interest  as  well. 

There  are  two  basic  formulations  of  the  DTOA  localization 
problem:  (a)  distance-differences  to  an  object,  such  as  the 
origin  of  a  plume  described  by  spatial  diffusions  [6],  are 


measured  from  known  locations,  and  the  problem  is  to  estimate 
the  location  of  the  object;  and  (b)  a  device,  such  as  a 
sensor  node,  receives  distance-differences  from  beacons  with 
known  locations,  and  the  problem  is  to  estimate  the  location 
of  the  sensor  node,  that  is  self-localization.  The  classical 
source  localization  problem  using  DTOA  measurements  has 
been  solved  using  two  general  approaches:  (i)  linear  alge¬ 
braic  solution  which  typically  involves  matrix  inversion  and 
solving  a  quadratic  equation  [2],  [7],  and  (ii)  intersection  of 
hyperbolic  curves  [8].  A  recent  overview  of  network-based 
localization  methods  may  be  found  in  [3]— [5],  [9].  In  general, 
the  quality  of  the  location  estimate  is  a  complex  function  of 
the  precision  with  which  the  underlying  numerical  operations 
are  implemented,  and  consequently,  there  is  no  apparent  and 
simple  way  of  relating  the  computations  to  the  “quality” 
of  location  estimate.  In  particular,  it  is  unclear  if  devoting 
more  computational  operations  would  increase  the  accuracy 
of  these  methods,  or  conversely  if  it  is  possible  to  reduce 
the  computations  without  drastically  affecting  the  quality  of 
location  estimate.  In  addition,  sensor  errors  can  have  drastic 
effects  on  DTOA  localization  methods.  For  example,  as  will 
be  shown  in  Section  IV  under  simple  random  noise  conditions, 
the  quadratic  equation  of  [2],  [7]  may  have  imaginary  roots 
thereby  rendering  the  method  incomplete.  Also,  numerical 
instabilities  may  arise  in  the  computations  implemented  with 
low  precision  operations  wherein  matrix  inversions  needed  for 
linear  algebraic  methods  may  become  ill-conditioned  resulting 
in  large  estimation  errors. 

The  underlying  geometric  nature  of  this  problem  has  been 
well-known  [1]  although  we  are  unaware  of  methods  that 
exploit  it  to  fine  tune  computations  as  done  in  several  com¬ 
putational  geometry  methods  [10],  [1 1] 1 .  We  present  a  com¬ 
putational  geometric  method  for  DTOA  localization  based  on 
a  binary  search  on  an  algebraic  curve  defined  by  a  distance- 
difference  function.  We  exploit  the  monotonicity  of  the  di¬ 
rectional  derivative  of  the  other  distance-difference  on  it  to 
support  the  binary  search.  The  computational  complexity  of 
this  method  is  0( log2  1/7),  where  the  computed  solution  is 

'The  term  triangulation  has  also  been  used  in  the  context  of  computational 
geometry  and  kinetic  data  structures  [12]  to  refer  to  the  decomposition  of 
planar  regions  into  triangles  which  is  quite  different  from  its  usage  in  this 
paper. 
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Fig.  1.  Regions  of  monitoring  area:  (a)  inside  region,  (b)  apex  region,  (c) 
bottom  region,  (d)  top  left  and  right,  (e)  bottom  left  and  right. 


guaranteed  to  be  within  a  7-precision  region  centered  at  the 
actual  solution;  by  a  minor  modification  to  the  termination 
condition  this  region  can  be  changed  to  [—7,7]  x  [— 7,71- 
box  centered  at  the  actual  solution.  Alternatively,  by  fixing 
the  number  of  operations  to  k,  one  can  achieve  the  precision 
7  =  0  ^2  .  This  method  is  robust  with  respect  to  distance 

measurement  errors:  (i)  7  is  of  the  same  order  of  magnitude 
as  errors  in  distance  measurements;  in  methods  that  involve 
division  operations  such  guarantees  cannot  be  made;  and  (ii)  it 
is  complete  in  that  it  will  always  return  an  answer,  even  under 
random  measurement  errors.  This  method  is  a  generalization 
of  the  DTOA  localization  method  in  [6]  proposed  as  a  part 
of  plume  identification  when  the  source  is  inside  the  acute 
triangle  formed  by  sensors.  In  our  case,  the  object  can  be 
located  anywhere  in  the  monitoring  region.  In  addition,  we 
also  provide  a  detailed  analysis  of  the  underlying  computation 
and  the  proof  of  the  required  monotonicity  property  of  the 
underlying  directional  derivative. 

This  paper  is  organized  as  follows.  We  describe  our  ge¬ 
ometric  DTOA  triangulation  method  in  Section  II.  We  prove 
the  correctness  of  the  method  by  establishing  the  monotonicity 
properties  of  the  underlying  directional  derivative  in  Section 
III.  We  present  simulation  results  in  Section  IV. 

II.  Geometric  DTOA  Method 

We  are  given  three  sensors  Si,  i  =  1,2,3  located  at  ( Xi,yi ), 
i  =  1,2,3.  For  any  point  P  =  (x,y)  in  the  plane,  we  define 

d(P,  Si)  =  y(x  -  x*)2  +  (2/  -  2/i)2  and 

A  ij(P)=d(P,Si)-d(P,Sj), 

for  i,  j  =  1,  2,  3.  We  consider  the  DTOA  localization  problem 
of  estimating  the  location  of  a  source  S  from  the  measurements 
of  Ai2(S)  and  A 13 (S'),  given  by  #12  and  #13,  respectively. 
As  we  move  P  from  Si  to  Sj  along  the  line  segment  SiSj, 
A ij(P)  varies  monotonically  and  linearly  from  —  d(Si,  S2)  to 


algorithm  geometric_DTOA(^i2  A13); 

begin 

L  (£12,2/12)  <—  intersection  point  of  L12  with  S1S2', 

2.  Xxi  set  of  X-coordinates  of  intersections  of  L12  with  Si  S3; 

3.  2x2  set  of  X-coordinates  of  intersections  of  L12  with  S2S3; 

4.  2x  «—  {Dxi,  £12, 2^x2}  LJ2xi  LJ2x2; 

5.  ls  <-  sort  (2xi ); 

6.  let  Is  =  {x(i),a;(2),...,x|Is|}  and  {yw,y(2),  ■  ■  ■ ,  Y\is\} 

7.  be  the  corresponding  Y -coordinates; 

8.  for  i  =  1, .. .  |XS|  -  1  do 

9.  ®(i)  <-max(i(i),flxi);  <-  min(x(i),Dx2); 

10.  y(i)  «-  max(t/(i),Dyi);  y(A)  <-  min (t/(i),  £>y2); 

11.  5^0; 

12.  for  i  =  1,. ..  |ZS|  -  1  do 

13.  □  <-  region_sign^a!(,>+^i+1?  ,^,r/(i+1)); 

14.  S  <-  S  1J{  locate.Li3  (x(i),a;(i+1),r/(i),s/(i+1),  □)}; 

15.  return  S; 

end 


d(Si,S2),  and  equals  0  at  the  bisector  point.  We  consider  the 
locus  of  points  defined  by 


LiJ(5)  =  {P\Aij(P)=5} 


which  is  described  by  the  algebraic  equation  d(P,  Si)  — 
d(P,Sj)=6. 

We  consider  the  generic  configuration  shown  in  Figure  1 
such  that  5 12  <  0  and  S 13  <0.  In  the  next  section,  we  show 
that  any  configuration  can  be  transformed  into  the  generic 
one  shown  in  Figure  1  that  consists  of  seven  regions.  In 
each  of  the  seven  regions  Ai3(.)  monotonically  varies  on 
L i2(-);  in  particular,  it  monotonically  decreases  inside  the 
apex  and  bottom  regions  and  monotonically  increases  in  all 
other  regions  as  will  be  proven  in  the  next  section.  The  overall 
strategy  to  estimate  S  =  (x,y)  is  to  perform  a  binary  search  on 
I/12 (•)  to  locate  S  =  (x,y)  such  that  |Ai3(5)  —  Sis\  <  7.  For 
any  point  P  in  the  plane,  let  Rp ;7  be  the  j -precision  region, 
which  corresponds  to  a  “distorted  box”  centered  at  P,  whose 
sides  are  formed  by  displaced  hyperbolic  curves  as  shown  in 
Figure  2.  The  above  condition  implies  that  S  E  Rs, 7. 

As  we  move  P  along  L^i^)  in  one  direction,  A  13(F) 
varies  monotonically  within  each  region.  The  basic  idea  is 
to  utilize  this  monotonicity  of  A 13  (P)  to  support  a  binary 
search:  repeatedly  compute  a  P  along  L  12(^12)  until  we 
reach  S  E  Rs, 7.  The  details  are  presented  in  algorithm 
geometric _DTOA(^i2,  £13). 


Fig.  2.  Rs ;7  is  a  “distorted”  box  corresponding  to  7-precision  region 
centered  at  S. 


algorithm  region_sign(x,  vltVr)\ 

begin 

1.  y  <-  locate_Li2(y/v,  y/?, a?); 

2.  (7i  «-  sign(y  -  yi  -  (y2  -  yi)0  -  xi)/(x2  -  xi)); 

3.  cr2  <-  sign(y); 

4.  (73  <—  sign(y  -  yix/xi); 

5.  if  ((7i,  (72,  (73)  =  (+,+,+)  or  then  return(<); 

6.  else  return(>); 

end 


Let  [Dx1,Dx2]  x  [LVi,Dy2]  be  the  monitoring  region 
within  which  S  is  to  be  localized.  The  basic  idea  of  algorithm 
geometric _DTO A  is  identify  individual  segments  of  L 12  that 
are  entirely  contained  in  single  regions  shown  in  Figure  1, 
and  perform  a  binary  search  on  L 12  with  L13  as  the  objective 
function  within  the  region.  The  correct  sign  (>  or  <)  for 
the  search  within  the  region  is  supplied  by  the  function 
region_sign(##,  2/l,  hr)  as  follows.  This  function  computes  a 
point  P  =  (x,y)  on  L12  (line  1)  and  evaluates  the  triple  of 
signs  by  substituting  x  and  y  into  the  equations  of  lines  S1S2, 
Si  S3,  and  S2S3  in  lines  2,  3  and  4  respectively.  It  then  returns 
<  if  the  computed  triple  matches  that  of  apex  or  bottom  region 
in  line  5,  and  returns  >  otherwise  in  line  6. 

In  algorithm  geometric_DTOA,  the  individual  segments  of 
L 12  that  are  entirely  contained  in  single  regions  of  Figure  1  are 
identified  in  lines  1-3,  and  are  arranged  in  ascending  order  of 
X-coordinates  of  the  end  points  of  these  segments  in  lines  5-7. 
The  segments  that  lie  outside  the  monitoring  region  are  con¬ 
strained  to  be  within  [Dx1 ,  Dx2\  x  [Dyi,Dy2\  in  lines  8-10. 
The  x  and  y  coordinates  of  points  on  L12  within  each  region 
are  bounded  within  the  intervals  [x^ ,  #(i+1)]  and  [y^ ,  y^^\, 
respectively  in  each  iteration  of  lines  13-14;  note  that  each  of 
corresponds  to  L12  intersecting  the  lines  through  two  of 
Si,  S2  and  S3,  or  the  end  points  of  [Dxi,Dx2\-  A  binary 
search  on  L 12  with  L13  as  the  objective  function  is  carried  out 
within  each  region  by  calling  locate _L  13  (.)  in  line  14  and  the 
returned  points  are  accumulated  in  S.  This  algorithm  returns 
S  which  contains  one  or  two  candidates  for  S.  For  simplicity, 
S  e  S  on  the  boundary  of  the  monitoring  region  is  interpreted 
as  either  as  a  source  lying  on  the  boundary  or  outside  the 
region.  Computation  of  (#12,2/12)  in  line  1,  Xxi  in  line  2  and 
Xx2  in  line  3  is  carried  out  by  a  binary  search  on  the  line  S1S2, 
Si  S3,  and  S2S3,  respectively  with  L12  as  objective  function. 


algorithm  locate_Li3(a;L,  xR,  yL,  yR,  □); 

begin 

1-  x  <-  (xL  +  xr)/2- 

2.  y  <—  locate_Li2(x,  yL,  Vr)\ 

3.  P  —  (x ,  y)? 

4.  if  |Ai3(P)  —  <5i3 1  <  7  then  return(P); 

5.  else  if  (Ai3(P)n<Si3)  then 

6.  locate-Li3(x,xR,yL,yR,  □); 

7.  else 

8.  locatQ.Li3(xL,x,yL,yR,  □); 

end 


As  shown  in  algorithm  locate  _L  13  (#l,  xr,  yL,  yR-,  LI), 
within  each  region  the  search  is  two-dimensional;  it  returns 
P  such  that  |Ai3(P)  —5 13 1  <7  that  is  contained  in  the  box 


algorithm  locate_Li2(x,  y^); 

begin 

i-  y  <-  (vl  +  vr)/2;  p  =  {x,y)\ 

2.  if  |Ai2(P)  —  S i2|  >  7  then  return(y); 

3.  else  if  Ai2(P)  <  <5i2  then  locate_Li2(x,  y,  y#); 

4.  else  locate_Li2(x,  y^,  y); 

end 


[xL,xR\  x  ['!iL,  yn\-  First,  in  lines  1-3,  a  point  P  =  (x,y)  on 
1/12  is  located  at  mid  #-value  computed  in  line  1  by  performing 
a  binary  search  for  y  by  algorithm  locate_Li2(#,  yL ,  Vr)  in  line 
2.  Once  A  13(F)  is  computed,  the  #-range  is  suitably  halved, 
and  this  process  is  recursively  carried  out  until  the  required 
precision  7  is  reached  in  lines  4-8. 

Complexity  of  computation  of  each  element  of  Xxi  and 
Xx2  is  0(log(l/7)),  and  each  has  at  most  2  points.  In 
computing  S  =  (x,y),  there  are  altogether  0(log(l/7))  calls 
to  locate _Li3(.),  and  each  call  invokes  locate_Li2(.X  which  in 
turn  has  a  complexity  of  0( log(l/7)).  Thus  the  complexity  of 
algorithm  geometric  DTOA(^i2,  73)  is  0(log2(l/7)),  which 
can  be  adapted  by  suitably  specifying  7.  We  have  \Is\  <  5, 
since  \Xxi  \  <  2  and  |Xx2 1  <  2,  and  thus  there  are  at 
most  four  invocations  of  locate _L  13.  If  the  number  of  basic 
computational  operations  are  fixed  k,  then  we  have  7  < 

O  (2-V^y 

III.  Monotonicity  of  Directional  Derivative 

In  this  section,  we  establish  the  correctness  of  the  method 
described  in  previous  section.  First,  any  given  configuration 
of  three  sensors  can  be  rotated  and  relabeled  such  that  vertex 
S 1  with  both  6 12  and  £13  are  negative  is  above  y  axis,  and  S3 
and  S2  are  aligned  along  #-axis;  note  that  S\  is  the  closest 
vertex  to  S  hence  always  exists.  Then  a  translation  ensures 
that  S 1  =  (#1,2/1),  S2  =  (#2,0)  and  S3  =  (0,0),  and  #1  > 
0;  yi  >  0;  #2  >  0  which  establishes  that  the  configuration  in 
Figure  1  is  generic. 

We  consider  five  separate  regions:  (a)  inside  triangle,  (b) 
top  apex,  and  (c)  bottom  region,  (d)  bottom  left,  and  (e)  top 
right  as  shown  in  Figure  1.  The  other  two  regions,  namely 
top  left  and  bottom  right,  are  “flipped”  versions  of  cases  in 
(e)  and  (d),  respectively,  and  can  be  similarly  handled.  We 
show  that  the  directional  derivative  of  Ai3(.)  along  the  curve 
I/12  (.)  is  monotone  in  each  of  these  regions:  it  is  positive  in 
regions  (a),  (d)  and  (e),  and  is  negative  in  regions  (b)  and  (c). 
The  correct  sign  of  the  region  for  the  search  is  supplied  by 
region_sign(.),  which  establishes  the  correctness  of  algorithm 
geometric  _DTOA. 

We  have  for  i  =  1,  2,  3, 

dd(S,Sj)  ( x-xj )  dd(S,Sj)  ( y-yi ) 

dx  d(S,Si )  an  dy  d(S,Si)' 

Then  directional  derivative  of  Ai3(P)  at  P  =  (x,  y)  on  the 
locus  .£<12(72)  =  {P|Ai2(P)  =  <>12},  for  any  72,  is  given 


by 


V  £l2(<5l2)Al3(-P) 


3Ai  2(P) 
dx 

gAi  2(P) 
dy 


(dA13{P) 
l  dec 


)  +(^) 


d(5,5i)  d(S,S3) 

y—yi _ y-y3 

d(S,S  i)  d(S,S3) 


K 


dA13(P) 

dx 

dAx  3(P) 
dy 


d(S,S  i) 
Ol 


d(S,S2) 

_ 2/~2/2 

d(5,5i)  d(S,S2) 


where  iT  = 
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v  s:y  , 


We  use  the  two 


following  basic  identities  extensively  in  our  derivations: 

1  —  cos  a  =  2  sin2  a/2 


cos  a  —  cos  (3  =  —  2  sin 


A.  Inside  Triangle 

In  this  case,  we  have  0  <  #  +  72  <  tt  and  #  >  71  as  shown 
in  Figure  3.  The  directional  derivative  is  given  by 
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We  have  0  <  71  +  72  <  7r,  0  <  #  +  72  <  tt  which  makes 
the  first  two  sin  terms  to  be  positive.  We  have  #  >  71  and 
0  <  #  <  7r.  Thus  — 7r/2  <  <  0,  which  makes  the 

third  cos  term  to  be  positive.  Hence  the  directional  derivative 
is  positive. 


B.  Top  Apex 

In  this  case,  we  have  0  <  #  +  71  <1 r  and  #  >  72  as  shown 
in  Figure  4.  The  directional  derivative  of  A  (Si,  S3)  on  the 
locus  {(#,  1/)  |  A  (Si,  S2)  =  <5i2 },  for  any  S12,  is  given  by 


Va(Si,S2)A(Si,  S3) 
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Fig.  4.  Source  S  =  (x,y)  is  located  in  the  apex  region. 
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We  have  0  <  71  +  72  <  7r,  which  makes  the  first  term  sin 
positive.  We  have  0  <  0  +  71  <7 r,  which  makes  the  second 
term  cos  positive.  Since  0  >  72  and  0  <  #  <  7r,  we  have 
— 7r/2  <  <  0’  which  makes  the  third  term  sin  negative. 

Hence  the  directional  derivative  is  negative. 


C.  Bottom,  Bottom  Left,  and  Top  Right  Regions 

(i)  Bottom  Region:  For  bottom  region,  the  derivation  is 
identical  to  the  case  of  apex  region:  the  conditions  0  < 
6  +  71  <7 r  and  0  >  72  are  valid  based  on  the  geometric 
conditions  specific  to  this  case  as  shown  in  Figure  5. 

(ii)  Bottom  Left:  The  case  of  bottom  left  is  identical  to  the 
bottom  region  except  that  0  +  71  >  7r  as  shown  in  Figure 
III-A,  which  makes  the  cos  term  negative,  and  hence  the 
directional  derivative  is  positive. 

(iii)  Top  Right:  The  case  of  top  right  region  shown  in  Figure 
6  is  identical  to  top  apex  region  except  that  2tt  >  #+71  > 
Q  +  a\  >  7 r,  which  makes  the  second  (cos)  term  negative, 
and  hence  the  directional  derivative  is  positive. 

Computational  results  indicating  the  signs  of  the  directional 
derivative  of  randomly  generated  sources  are  shown  in  Figure 
7. 


IV.  Simulation  Results 

We  simulated  a  network  of  three  sensors  on  a  [0, 100000]  x 
[0, 100000]  grid  such  that  S3  and  S2  are  located  at  (0, 0) 
and  (100000,0)  respectively.  Location  of  S%  is  randomly 
generated  on  the  line  segment  between  (0, 100000)  and 
(100000,100000).  Each  sensor  measurement  corresponds  to 
(1  +  f)r  where  r  is  the  actual  distance  from  sensor  to  source, 
and  /  is  uniformly  randomly  generated  in  the  interval  [0,  F] 
for  a  fixed  multiplicative  factor  F.  While  /  values  are  gen¬ 
erated  independently,  sensor  error  magnitude  is  proportional 
to  the  distance  from  the  sensor  to  the  source.  The  sensor 
errors  are  correlated  due  to  the  spatial  relationships  between 
the  sensor  locations;  a  source  close  to  S3  generates  a  small 
error  there  and  larger  errors  at  both  Si  and  S2,  which  are 


Fig.  6.  Source  S  =  (x,y)  is  located  in  the  top  right  region. 


Fig.  7.  Source  S  =  ( x,y )  is  randomly  selected,  and  the  sign  of 
Va(Si,S2)  A(Si,  S3)  is  computed. 


F 

imaginary  roots 
percentage 

0.01 

0.01 

0.02 

0.047 

0.03 

0.13 

0.05 

0.32 

0.10 

0.963 

TABLE  I 

Listing  of  percentage  of  imaginary  solutions  to  the  DTOA 

QUADRATIC  EQUATION  AS  A  FUNCTION  OF  THE  MULTIPLICATIVE  FACTOR 

F. 


located  farther  away.  From  these  measurements,  we  computed 
distance-differences  and  tested  DTOA  localization  methods. 

We  implemented  a  linear  algebra  based  method  of  [2],  [7] 
which  required  a  solution  to  a  quadratic  equation.  When  sensor 
errors  are  zero  (F  =  0)  this  method  accurately  estimated  the 
source  location.  However,  when  F  >  0,  this  method  became 
incomplete  as  the  quadratic  equation  had  imaginary  roots  for 
a  small  percentage  of  sources  as  shown  in  Table  1  based 
on  simulation  of  100,000  sources.  The  generated  sources  are 
shown  in  Figure  8(a)  which  are  uniformly  distributed  across 
[0, 100000]  x  [0, 100000]  grid.  For  the  case  F  =  0.05  in  Table 
1,  about  0.32%  of  the  sources  yielded  imaginary  solutions  to 
the  quadratic  equation,  and  these  sources  themselves  are  are 
concentrated  around  the  sensors.  On  a  related  note,  the  method 
of  [13]  accounts  for  random  errors  that  are  independent 
Gaussian,  and  hence  is  not  directly  applicable  to  this  case. 

Results  of  our  method  are  shown  in  Figure  9  for  10 
different  values  of  7  with  7min  =  7.856613  and  7max  = 
29.011267.  This  method  is  complete  in  that  always  returned 
the  precision  region.  When  sensor  errors  are  zero,  this  region 
always  included  the  source.  If  sensor  errors  are  non-zero, 
this  region  did  not  always  include  the  source  depending  on 
the  value  of  7.  For  each  value  of  7,  we  list  fs  which  is 
the  fraction  of  the  sources  that  were  not  included  in  the 
precision  region.  As  expected  smaller  value  of  7  resulted  in 
more  missed  sources  and  a  large  enough  value  always  included 
the  source.  Also,  for  each  value  of  7,  we  are  list  the  maximum 
number  of  steps  needed  in  computing  S  in  Figure  9  over  the 
100,000  sources.  This  method  is  implemented  in  C  on  Linux 
workstation  (Opteron  processor),  and  the  typical  execution 
times  of  algorithm  geometric _DTO A  for  datasets  are  under 
a  second. 


V.  Conclusions 

We  presented  a  computational  geometric  method  for  the 
problem  of  triangulation  in  plane  using  measurements  of 
distance-differences.  This  problem  has  been  extensively  stud¬ 
ied  in  the  past  and  several  solutions  have  been  deployed,  and 
our  re-examination  is  motivated  in  part  by  the  requirements 
of  low  power  sensor  nodes.  Our  method  is  computationally 
efficient  and  adaptive  as  well  as  robust  with  respect  to  mea¬ 
surement  and  computational  errors.  This  method  is  particularly 
suited  for  deployment  in  nodes  that  adapt  their  computations 
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(b)  sources  that  yielded  imaginary  roots 

Fig.  8.  Sources  that  yielded  imaginary  roots  in  DTOA  localization  methods 
based  on  quadratic  equations  are  concentrated  around  the  sensors. 
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Fig.  9.  Fraction  of  sources  lying  outside  Rg  corresponding  to  different  7 
values. 
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in  response  to  power  budgets.  This  method  can  also  be  applied 
when  distance  measurements  are  available,  and  can  offer  simi¬ 
lar  advantages  over  the  linear  algebraic  methods  that  are  often 
used  for  triangulation  based  on  distances.  Furthermore,  by 
computing  distance-differences  from  distance  measurements, 
this  method  would  be  less  susceptible  to  one-sided  bias  errors 
in  distance  measurements.  This  is  particularly  useful  in  certain 
self-localization  tasks,  where  a  single  sensor  is  employed  to 
measure  distances  to  reference  beacons. 

This  paper  is  only  a  step  towards  utilizing  computational 
geometric  methods  for  solving  localization  problems.  It  would 
be  of  future  interest  to  consider  extensions  of  this  method 
for  cases  where  more  than  three  sensors  are  deployed  and 
multiple  measurement  sets  are  provided.  It  would  also  be 
interesting  to  see  if  the  proposed  method  can  be  extended 
under  random  noise  models.  Also,  multiple  path  effects  are 
not  considered  in  this  paper,  and  it  would  be  of  interest  to 
explore  such  extensions.  For  the  special  case  when  Si,  S2  and 
S3  form  an  acute  triangle,  a  training  method  was  proposed  in 
[6]  wherein  the  localization  method  can  be  trained  in- situ  to 
account  for  sensor  correlations.  The  current  method  can  be 
similarly  employed  but  the  training  procedure  is  likely  to  be 
more  involved.  It  would  be  of  future  interest  to  explore  the 
“tracking”  ability  of  this  method  by  repeatedly  executing  it  on 
a  stream  of  distance-difference  measurements  corresponding 
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